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Abstract.  For  a  scalar  conservation  law  Q  =  with  of  constant  sign,  the  first  order 

upwind  dilTerence  scheme  is  a  special  case  of  Godunovas  rnotliod.  The  method  is  equivalent  to 
solving  a  sequence  ofRiemann  problems  at  each  stop  and  averaging  the  resulting  solution  over  each 
cell  in  order  to  obtain  the  numerical  solution  at  the  next  time  level.  The  diiTcrcnce  scheme  is  stable 
(and  the  solutions  to  the  associated  sequence  of  Riemann  problems  do  not  interact)  provided  the 
Couraiit  number  u  is  less  than  1.  Ry  allowing  and  explicitly  handling  such  Interactions,  it  is  possible 
to  obtain  a  generalized  method  which  is  stable  for  u  much  larger  than  1.  fn  many  cases  the  resulting 
solution  is  considerably  more  accurate  than  solutions  obtaine<l  by  other  numerical  metho<ls.  In 
particular,  shocks  can  be  correctly  computed  with  virtually  no  smearing.  The  generalized  method 
is  rather  unorthodox  and  still  has  some  problems  associated  with  it.  Monctholcss,  preliminary 
results  arc  quite  encouraging. 
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1»  Introdiiction. 


A  scalar  conservation  law  Is  a  partial  difTcrential  equation  of  the  form 


Ut(x,t)  =  /(u{i,t))* 


(l.l) 


where  «  :  IR  X  [0, oo)  — *  IR  and  /  ;  IR  — >  HI.  When  u  is  a  density  and  /  a  Ilux,  (1.1)  states  that 
the  integral  of  u  over  some  interval  [xx,x^z)  changes  only  due  to  the  flux  through  the  endpoints, 

/*2  Q 

^(/(w(a:) 

=  f{u.(x2,t))  -  f(u(xi,t)). 

The  theory  of  conservation  laws  is  described,  for  example,  by  Lax(7l  and  Whithatn[10]. 

The  most  basic  problem  for  a  conservation  law  is  the  /^iema/m  prob/em,  which  is  (l.l)  together 
witli  piecewise  constant  initial  conditions  with  a  single  discontinuity, 


^  r 

dt  7x1 


X  <  Zu 
X  >  Xq. 


(1.2) 


The  solution  of  a  Ricmann  problem  can  often  be  found  analytically.  Various  numerical  schemes 
for  solving  (l.l)  with  arbitrary  initial  data  arc  based  on  solving  a  sequence  of  Ricmann  problems 
exactly.  As  usual  with  linite  difference  schemes,  we  choose  a  spatial  stepsize  h  and  a  tirnostep  k 
and  set  xy  =  jh^  tn  =  nk.  Wc  then  approximate  ti(xy,in)  by  u” .  But  now  wc  view  the  discrete 
solution  }  not  as  an  approximation  to  some  smooth  function  but  rather  as  a  representation  of 
a  step  function 

u{x,t„]  =  u”  for  a:y_,/2  <  a:  <  a:,  +  i/2  (1-3) 

where  Each  discontinuity  defines  a  single  Riernann  problem.  If  the  time-step 

k  is  sulliciently  small,  the  solutions  to  these  Riernann  problems  do  not  interact  in  time  k  and 
the  solution  to  (1.1)  with  initial  conditions  (1.3)  can  be  computed  exactly.  This  is  the  case  if  the 
Cotinint  nftmbcrf 

k 

!/  =  T  wax  /'(v(x,  tn)), 
n  X 

is  less  than  1.  In  order  to  continue  this  process  wc  must  project  the  exact  solutions  of  the  associated 
Ricmann  problems  onto  the  computational  grid.  This  mapping  can  be  defined  in  various  ways. 
Setting  to  the  average  value  of  u(x,<n+i)  over  the  cell  (x,_i/2,  Xy^i/2)» 

=  h  /  u(x,i„+i)dx 


gives  Godunov’s  met/iod(5].  Setting  to  the  value  of  u(x,in+i)  ‘''t  some  randofnly  chosen 

point  in  (xy_i/2» '•*^j4*i/2)  gives  Chorin’s  version  of  Gfimm's  sebome,  also  known  as  the  random 
choice  metbod(lI(3].  This  method  has  the  advantage  that  shocks  and  contact  discontinuities  in  the 
solution  remain  sharp.  However,  their  positions  arc  no  longer  correct  in  general,  although  they 
arc  in  an  average  sense.  Both  of  these  methods  can  also  be  used  more  generally  for  systems  of 
conservation  laws. 

The  aim  of  the  present  paper  is  to  show  how,  for  scalar  conservation  laws,  the  restriction 
1/  <  1  can  frequently  be  dropped  by  explicitly  handling  the  interactions  of  the  solutions  to  the 
various  Riernann  problems.  This  has  the  obvious  a<l vantage  of  rcriuiring  less  work  to  advance  the 


solution  to  a  given  lime.  It  also  turns  out  to  bo  more  accurate  in  many  situations.  In  particular, 
Tor  large  values  of  t/,  sharp  shocks  can  be  maintained  even  when  using  (]o<Iunov-type  averaging. 
The  problems  associated  with  random  choice  methods  can  thus  be  avoided.  The  resulting  scheme 
is  conservative  and  shock  locations  arc  correct,  at  least  Tor  problems  i/jvolving  only  shocks.  The 
procedure  described  below  has  also  given  very  good  results  on  numerous  problems  involving  shock- 
rarefaction  interactions,  but  there  are  still  some  diflicultics  to  be  resolved  in  this  area.  These  will 
he  discussed  in  section  5. 

In  some  ways  the  method  is  best  thought  of  as  a  generalization  of  the  first  order  upwind 
dilfercnce  scheme  for  a  constant  cocificicnt  linear  problem.  Recently  Roc[9]  has  advocated  viewing 
difference  schemes  for  such  equations  in  an  “increment  form"  in  order  to  generalize  them  to  new 
methods  for  systems  of  conservation  laws.  We  will  begin  in  the  same  framework,  and  generalize  in 
a  dilTerent  direction. 

2.  The  constant  coefficient  linear  problem. 

We  begin  by  considering  the  constant  coelTicient  linear  advection  equation ' 

tij  =  cti*  c  >  0.  (2.1) 


We  define 


u  =  ckjh,  the  Courant  number. 

Any  explicit  linear  two-levcl  difference  scheme  can  be  written  in  increment  form  as 


(2.2) 


for  some  choice  of  the  coefllcicnts  7,-.  The  scheme  is  at  least  first  order  accurate  if 


Further  conditions  on  the  7^  give  higher  order  methods.  Sec  Roc[9]  for  a  more  general  discus.sion. 
We  now  think  of  implementing  (2.2)  in  the  following  way:  for  each  j  wo  compute  i/Ay  ,  split  it  up  into 
fractions  proportional  to  7^,  and  add  the  ith  fraction  to  Thus  rather  than  interpreting  (2.2) 

as  a  formula  for  computing  we  can  take  it  as  a  prescription  for  distributing  the  incremental 

data  over  nearby  meshpoints. 

The  choice  of  coefTicienls  7^  can  be  motivated  in  the  following  way.  For  simplicity  consider  a 
single  Riemann  problem 


j  <  J 
j>J. 


(2.3) 


This  is  interpreted  as  a  discretization  of 


X  <  Xj+x/i 

X  >  Xj  +  x/2- 


At  time  tn+i  the  solution  of  this  Ricrnann  problem  is 


X  <  XjH-l/2  -  ''A 
^  ^  ^J+\/2  - 
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We  would  like  our  approximation  to  be  some  discrete  representation  of  this  jump.  If  i/  =  1, 
the  jump  is  at  xy_x/2  ^nd  is  represented  by 


j  ^  J  (u  — 


(2.4) 


On  the  other  hand,  if  <  1,  the  true  position  of  the  shock  can  be  represented  by  smearing  the 
jump  over  more  meshpoints.  One  possiblility  is 


j<J 

+t/{P  —  q)  3  —  J 

j>J. 


(2.5) 


This  agrees  with  (2.4)  for  t/  =  1  and  has  been  chosen  so  that  the  discrete  conservation  law  holds, 


Note  that  each  Uj  is  the  average  value  of  u(x,  t)  over  the  interval  (ij_i/2)  Xj-|_i/2l* 

For  arbitrary  {u"}  we  apply  (2.5)  to  each  discontinuity  in  the  step  function  (1.3).  This  yields 
the  standard  hrst  order  upwind  difference  scheme 

=  u,"  (2.6) 


obtained  from  (2.2)  by  taking  7o  =  1  and  all  other  7*  —  0. 

Other  conservative  representations  of  a  jump  at  X/-1-1/2  —  are  also  possible.  These  lead  to 
other  special  cases  of  (2.2).  For  example,  replacing  (2.5)  by 


u 


; 


'a 

a+ii/(l  +iy){l3  —  a) 
0 +lu{l  —  i/)(P  —  a) 

y0 


j  <j-i 
3  =  J 
J  =  J  +1 

3>J +2 


(2.7) 


leads  to  the  general  scheme 


+^‘'(1  +‘^)A? 


which  is  the  Lax-Wendroff  scheme  for  (2.1).  For  smooth  solutions  this  is  often  preferable  to  (2.6), 
being  second  order  accurate.  For  a  single  time-step  on  the  Riemann  problem,  however,  (2.7)  is 
clearly  inferior  to  (2.5)  since  the  discontinuity  has  been  smeared  over  more  mesh  points  than 
necessary  and  since  the  approximation  is  no  longer  monotonic. 

We  will  call  (2.5)  the  optimal  representation  of  a  jump  at  x/^-i/2  —  since  it  is  monotonic 
and  involves  as  few  mesh  points  as  possible. 

Implementing  (2.6)  can  now  be  described  as  follows:  interpret  the  data  {u”}  as  being  a  sequence 
of  Riemann  problems,  solve  the  Riemann  problems  exactly,  and  then  represent  the  resulting 
discontinuities  on  the  finite  mesh  via  (2.5).  This  viewpoint  will  prove  fruitful  when  dealing  with 
more  general  conservation  taws. 
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S.  Upwind  differencing  with  </  >  1. 

Adopting  the  viewpoint  proposed  in  the  previous  section  allows  us  to  extend  the  scheme  (2.6) 
for  use  with  large  time  steps.  Attempting  to  apply  (2.6)  directly  with  2/  >  1  is  unstable.  Instead 
we  again  interpret  the  discrete  data  as  a  representation  of  a  step  function  and  solve  a  sequence  of 
Riemann  problems.  Again  consider  the  single  discontinuity 


j<J 

j  >  J. 


The  discontinuity  A j  =  /?  —  a  at  xj^i/2  propagates  under  (2.1)  to  Xj^i/2  —  uh  at  time  fn+i- 
Let  fjL  =  [t/\,  the  int^er  part  of  2/.  Then  x /4.1/2  ^  x The  solution  to 
the  Riemann  problem  is 


X  <  Iy_^4.i/2  —  (1/  —  n)h 

X  >  ly-M+i/J  —  (1/  —  ti)h. 


(3.1) 


which  we  must  now  represent  on  the  discrete  mesh.  By  again  averaging  the  solution  over  each  cell 
we  obtain  the  representation 


fa  i  <  /  -  M 

=  <  a  +(^  —  —  a)  j  =  J  —  n  (3.2) 

U  j  >  J  —  fX. 

This  is  also  the  natural  generalization  of  (2.6)  if  we  read  (2.6)  not  as  **add  t/  times  to  the  next 
mesh  point  over’’  but  rather  as  “add  Ay  to  the  next  u  mesh  points  over”.  Applying  (3.2)  at  each 
discontinuity  gives  the  following  scheme  for  general  {«y  }■ 

=  u;  +  a;  +A7+,  +• . .  +a;^^  +(1.  -  m)a"+,+i. 

Most  of  these  terms  cancel,  leaving 

=  (1  _  (^  _  m))u^+^+1  +{<'  - 

This  is  simply  the  method  of  characteristics  with  linear  interpolation. 

4.  General  scalar  conservation  laws. 

Now  consider  a  scalar  conservation  law 


“t  =  /(u),.  (4.1) 

We  will  always  assume  that  /''(u)  has  constant  sign  for  all  u  of  interest.  Suppose,  to  begin  with, 
that  >  0  and  that  u(x,0)  is  nondecreasing. 

Again  let  A”  =  —  Uy  and  set 

c7  =  (/(«?+,)-/(«?))/ a;  for 
u;  =  c;k/h. 
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Each  Cy  is  the  propagation  speed  of  the  discontinuity  Ay  according  to  the  Rankine-Hugoniot 
relation.  We  do  not  need  to  define  if  Ay  =  0;  there  is  no  need  to  propagate  jumps  of  height 
zero.  In  general  the  discontinuity  Ay  is  now  a  shock  which  travels  at  speed  Cy  and  should  propagate 
grid  points  in  one  time  step.  For  concreteness  assume  that  it  is  propagating  to  the  left.  As  a 
first  attempt  at  generalizing  the  procedure  of  section  3  we  try  adding  A]^  to 

and  (i/y  —  /iy)A”  to  Uy_^n.  This  gives  a  conservative  scheme,  since 

>  j  > 

so 

E  -  E  K  =  *  -  /(“")) 

=  k(f{+Qo)—  /(— oo)). 

Again  this  scheme  can  be  interpreted  as  the  method  of  characteristics  with  linear  interpolation. 
The  problem  is  that  when  some  i/y  >  1  we  may  be  allowing  characteristics  to  cross,  i.e.  some 
shocks  may  propagate  right  past  their  neighbors,  rather  than  coalescing  Into  a  single  shock  as 
should  happen.  Aji  example  will  illustrate  this. 

Examp/e  4.1,  Consider  f{u)  =  Jw*  in  (4.1),  i.e.  Ut  =  uu*-  Take  the  initial  conditions 

fO  j  <J 

«»=  2  j  =  J  (4.2) 

U  j>  J 

and  k  =  2h.  Then 

Aj — 1  =  2  I  —  1  1  “  2 

=  1  4  =  5/2  1/^  =  5 

We  thus  add  A®r_,  to  the  two  mesh  points  u®/  ,  and  and  A®r  to  the  five  mesh  points 

u®,...,u®^,.w4  find  that 

fO  i  <  J  -  4 

uj  =r  h  <  j  <  J«3  (4,3) 

U  j>J—2 

This  is  illustrated  in  Figure  4.1a.  Rather  than  coalescing  into  a  single  shock,  A®  has  “passed  right 
by^  The  resulting  is  not  the  optimal  representation  of  the  true  solution  at  time  ti, 

although  it  is  a  conservative  representation.  (The  true  solution  is  a  single  shock  at  ^ 

Remarkably,  the  shock  does  not  become  further  smeared  out  at  subsequent  times.  In  fact, 
starting  from  (4.3)  we  compute  that 

f 0  j<J-7 

u*  =  h  i  =  j-7 

U  j>  J~7 

as  shown  in  Figure  4.2b.  In  this  case  the  shock  has  sharpened  and  is,  in  fact,  the  optimal 
representation  of  the  true  solution:  a  shock  of  height  3  at  xj^s  ^  It  has  precisely  the  same 
shape  as  u®  and  at  future  times  u®,u^,u^, . . .  will  alternate  between  the  shape  of  and  that  of 
u®,  always  being  a  conservative  and  monotonic  representation  of  the  true  solution  (although  not 
always  the  optimal  representation).  In  this  example  the  shock  is  never  badly  smeared  and  the 
solution  is  quite  acceptable  for  all  tn*  If  we  had  taken  k/h  much  larger,  however,  the  resulting 
solution  would  have  been  badly  smeared,  at  least  for  some  tn* 
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This  problem  can  be  greatly  alleviated  by  explicitly  handling  the  coalescing  of  shocks.  This 
is  not  as  difficult  as  it  sounds  due  to  the  following  result. 

Proposition  4.1.  Consider  the  conservation  /aw  ut  =  (/(u)),  with  /"(u)  >  0  for  a  <  u  <  7 
and  initial  conditions 

fa  X  <ii 

u{x,o]  =  \p  ei<^<6 

h  I  >  ^2 

for  some  $1  <  ^2#  a  <  ^  <  7.  Let  tc  be  the  time  at  which  the  two  shocks  coalesce.  Consider  the 
same  problem  Vt  =  (/(v))*  with  initial  conditions 


X  <  (o 

>  Co 


where  $0  =  ((fi  --  a)$i  +(7  —  —  a).  Then,  for  t  >  U, 


u(i,  t)  =  v(x,  t)  Vx. 


The  same  result  holds  if  a  >  0  >  'y  and  f"(u)  <  0  for  7  <  u  <  at. 

This  result  is  an  immediate  consequence  of  conservation.  The  implication  is  that  whenever  two 
adjacent  shocks  and  A}  are  going  to  coalesce  at  some  time  between  tn  and  (i.e.  when 

kc'}  >  kc^_j  +h)  we  can  simply  merge  them  into  a  single  shock  A  =  A^__j  +  Ay  at  time  tn  and 
propagate  this  merged  shock.  At  time  tn»  A  is  located  at  x  =  (A”_jXy_i/2  +AyX y^i/2)/A  and 
it  propagates  at  speed  c  =  —  /(^y— i)-  l^ben  optimally  represent 

the  resulting  shock  solution  at  time  tn+i  on  the  discrete  grid.  Applying  this  procedure  to  example 
4.1  gives 

{0  j  <  J  —  3n 

2  y  =  y  —  3n  Vn 

3  j  >  /  —  3n 

which  is  always  the  optimal  representation  of  the  true  solution. 

For  general  data  {u” }  it  may  be  necessary  to  merge  several  shocks  at  a  time.  Since  the 
propagation  speeds  of  the  shocks  change  as  they  coalesce,  it  is  not  always  easy  to  write  out  the 
appropriate  equivalent  problem  a  priori.  The  following  procedure  can  be  used,  assuming  A”  is 
nonzero  for  a  finite  number  N  of  points  Xj.  Let  T  =5  {T*}  be  the  set  of  ordered  triples 


Ti  =  (xi,u^  ,u+) 


where  each  x,  is  the  location  of  a  shock,  a  discontinuous  jump  in  u  from  the  value  u~  on  the  left 
to  on  the  right.  The  T*  are  ordered  from  left  to  right  so  that  x,  <  x*,^!.  For  each  triple  we 
define 

ft  =  (/(«,'^)  -  /(“.“))/{“.■’■  -  »*r) 

A,  =  u+  -  u- . 

Initially  the  T*  are  chosen  as  the  triples  corresponding  to  nonzero  AJ*.  We  wish 

to  replace  T  by  a  new  set  of  shocks  5  which  do  not  interact  in  time  k  and  which  yield  the  same 
exact  solution  at  time  ^.he  original  set.  This  can  be  accomplished  by  Algorithm  4.1,  which 

works  from  left  to  right  and  employs  a  stack  of  triples  St  for  *  =  1,2,.. .,  /.  The  rightmost  shock 
on  the  stack  is  represented  by  S/,  which  is  on  top  of  the  stack.  The  stack  has  the  property  that 
at  the  beginning  of  each  iteration  the  shocks  represented  by  the  5^  would  not  interact  in  time  k  if 
there  were  no  disturbance  to  the  right  of  X/  (i.e.  if  u  =  for  x  >  x/).  Each  new  shock  is  put 
on  top  of  the  stack  and  then  is  merged  as  necessary  with  shocks  already  on  the  stack. 
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ALGORITHM  4.1. 

/:=  1 
Si  :=  Ti 

for  j  :=  2,3,...,  N  do 
/:=/+! 

Si  :=  Tj 

while  0  <  (i/  —  X/_i)/(6/  —  c/_i)  <  A;  do 

remove  5/  and  5/— i  from  the  stack  and  replace  them  by  a  single  triple  with 

X  =  +A/) 

u-  =  ujl.1 
u+  =uf' 

/:=/-! 

It  looks  as  if  this  algorithm  may  require  0(N^)  steps  in  the  worst  case,  because  of  the  nested 
inner  loop.  It  is  always  linear,  however,  because  the  innermost  step  merges  two  shocks.  Since  we 
begin  with  only  N  shocks  this  will  be  executed  at  most  N  times. 

Proposition  4.2.  For  any  sfcepsise  k,  Algorithm  4J  computes  the  exact  solution  to  (4.1)  for 
piecewise  constant  initial  data  with  a  finite  number  of  shocks. 

Proof.  Two  shocks  which  appear  to  interact  in  time  k  when  considered  in  isolation  must  also 
interact  in  time  k  in  the  presence  of  other  shocks.  Since  the  merging  procedure  is  associative,  the 
order  in  which  we  merge  them  is  irrelevant.  | 

In  general,  errors  will  be  introduced  interpreting  and  representing  the  solution.  For  "smooth** 
portions  of  the  solution  the  method  can  again  be  viewed  as  the  method  of  characteristics  with 
linear  interpolation.  The  error  in  a  single  time  step  is  thus  O(h^),  giving  0{h}/k)  global  error. 
A  true  discontinuity  in  the  solution  is  propagated  to  exactly  the  correct  location,  providing  the 
solution  is  properly  interpreted. 


5.  Rarefaction  waves. 


Now  consider  a  jump  for  which  Ay/"(u)  <  0.  Such  a  discontinuity  should  spread  out  into 
a  rarefaction  wave.  Unfortunately,  applying  the  procedure  of  section  4  directly  will  cause  A”  to 
propagate  as  a  sharp  discontinuity. 

In  dealing  with  shocks  we  sometimes  found  it  necessary  to  merge  several  discontinuities  into 
a  single  shock.  Rarefaction  waves  can  be  approximated  by  taking  the  opposite  approach.  The 
discontinuity  A^  is  broken  up  into  several  (say  m)  discontinuities  each  of  height  A”/m  and  all 
located  at  These  discontinuities  will  then  travel  at  different  speeds,  spreading  out  to 

approximate  the  true  rarefaction.  We  should  choose  m  large  enough  that  the  resulting  solution 
looks  smooth.  Since  the  true  rarefaction  spreads  over  meshpoints,  it  is 

reasonable  to  choose  m  as  an  approximation  to  this  quantity,  e.g.  as  the  integer  part  of  some  finite 
difference  approximation.  In  general  rarefaction  and  shock  waves  may  interact.  It  will  again  be 
necessary  to  handle  this  interaction  explicitly  in  order  to  obtain  good  results  for  large  values  of 
k/h.  These  interactions  can  also  be  handled  by  Algorithm  4.1,  at  least  approximately,  if  we  define 
m  triples  Ti,  t  s  1, 2, . .  .,m  corresponding  to  each  partitioned  discontinuity  with 


^  m  ^ 
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See  Figure  5.1  for  an  example  of  this. 


Example  5.1.  In  order  to  better  understand  this  shock-rarefaction  interaction,  consider  the  problem 
ut  =  with  initial  data 


3  =  J 
j  ^  */• 


(5.1) 


This  is  a  representation  of 


Xj^il2  <  X  <  X/4-1/2 

otherwise. 


(5.2) 


Take  k  ^  h.  Then  the  true  solution  at  time  k  is 

u(x,k)  =  Xy^.i/2  —  y/Sh  <  X  <  Xy^i/2  (5  3J 

lo  otherwise. 

The  rarefaction  has  overtaken  the  shock,  slowing  it  down  and  decreasing  its  strength.  In  the  spirit 
of  Proposition  4.1  we  can  attempt  to  find  a  new  set  of  initial  conditions  which  give  the  same 
solution  at  time  A;  without  any  interaction.  The  desired  initial  conditions  arc 

v{x,0)  =  (^  — V2h  <  X  <  Xj^i/2  (5.4) 

lo  otherwise 


This  pulse  is  wider  and  shorter  than  the  original  pulse  (see  Figure  5.1).  In  the  solution 
there  is  no  interaction  before  time  k  at  which  point  the  rarefaction  wave  just  catches  up  to  the 
shock.  For  t  >  k,  v(x,  t)  =  u(i,  t)  for  all  x.  Computing  the  exact  equivalent  set  of  noninteracting 
discontinuities  (such  as  (5.4))  for  a  general  scalar  conservation  law  with  step  function  initial 
conditions  would  be  difficult  but  could  be  done.  Taking  this  approach  and  then  solving  the  resulting 
Riemann  problems  exactly  would  yield  the  exact  solution  to  such  a  problem,  just  as  Algorithm  4.1 
gives  the  exact  solution  for  a  pure  shock  problem. 

Instead  of  following  this  course,  we  will  concentrate  on  the  more  easily  obtainable  approximate 
solution  given  by  Algorithm  4.1  together  with  the  previously  described  rarefaction  partitioning 
procedure.  Taking  m  =  4  in  the  rarefaction  gives  us  5  triples  Tt  in  the  set  T.  These  are 


Ti  =  (Xy— 1/2,  0,  4), 

72  =  (Xy-f.i/2,  4,  3), 

Ta  =  (xy+i/2,  3,  2), 

^4  =  (xy+i/2,  2,  1), 

^5  =  (Xy4.i/2,  1,  0). 

In  applying  Algorithm  4.1,  we  see  that  Ti  and  7*2  interact.  This  is  the  only  interaction  and  the 
resulting  set  S  consists  of  4  triples.  The  first  of  these,  is  derived  from  Ti  and  7*2  and  has 


Xi  =  ^4Xy_i/2  —‘Xj^i/2) 
=  xy  —  J/l 
up  =  0 
=  3. 
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r 


Figure  5.1  shows  the  initial  data  both  before  and  after  merging.  The  original  pulse  has  been 
replaced  by  a  wider  and  shorter  pulse  which  approximates  (5.4).  We  advance  the  solution  to  time 
k  by  moving  each  Sx  at  the  appropriate  speed.  The  final  positions  are 


Si  ■■ 

-  \h)  -  \h 

II 

H 

1 

1 

5j: 

I/+1/2  — 

=  2^/— 3 

53: 

^y+1/2  — 

= 

54: 

iy+1/2  — 

=  X/ 

This  is  shown  in  Figure  5.2a.  These  discontinuities  are  then  represented  on  the  mesh  in  the  usual 
way;  giving  the  solution 

12  i  =  7  -  2 

3/2  j=J-l 

1/2  j  =  J 

0  otherwise 

which  is  shown  together  with  the  exact  solution  in  Figure  5.2b. 

In  practice  it  is  necessary  to  modify  this  procedure  slightly  in  order  to  maintain  smoothness 
in  the  rarefaction  waves  over  several  time  steps.  Rather  than  putting  each  jump  in  the  partitioned 
rarefaction  at  we  spread  them  out  between  and  Xj-^i  by  setting 

=  *  =  1,2 . m. 

m  +1 

This  procedure,  together  with  Algorithm  4.1,  has  given  quite  satisfactory  results  for  a  variety 
of  problems.  It  must  be  stressed,  however,  that  Algorithm  4.1  may  not  handle  interactions  properly 
when  several  shocks  and  rarefactions  are  present.  Consider,  for  example  a  problem  involving  two 
shocks  followed  by  a  rarefaction.  It  may  be  that  the  two  shocks,  considered  in  isolation,  would 
coalesce  in  time  A:,  whereas  in  reality  the  the  rarefaction  catches  up  to  the  second  shock  and  slows 
it  down  sufficiently  that  no  further  interaction  occurs.  Algorithm  4.1,  which  works  from  left  to 
right,  would  produce  an  incorrect  result  which  may  be  quite  bad  for  large  Courant  numbers. 

Even  more  disturbing  is  the  fact  that  monotonicity  of  the  Xx  may  be  lost  when  merging  a 
shock  with  a  rarefaction.  If  A/  and  A/_i  have  opposite  signs  then  the  new  x  defined  in  Algorithm 
4.1  does  not  lie  between  X[  and  and  hence  may  lie  to  the  left  of  z/_2  or  to  the  right  of 

the  next  Currently  this  is  handled  by  aborting  execution  whenever  x/  <  X[^\,  This  has 

never  happened  for  moderate  values  of  the  Courant  number  on  “reasonable”  initial  data. 

6.  Numerical  results  for  ID  scalar  conservation  laws 

Before  proceeding  to  briefly  discuss  several  space  dimensions  and  systems  of  conservation  laws, 
we  pause  to  present  some  computational  examples. 

Figure  6.1  shows  several  sets  of  initial  data.  The  remaining  figures  show  the  results  of  applying 
the  method  to  problems  with  those  initial  conditions.  In  Figures  6.2  through  6.4,  the  true  solutions 
are  known  and  are  shown  as  solid  lines.  The  numerical  solutions  are  represented  by  the  squares. 

The  other  examples  begin  with  more  complicated  initial  data  and  use  periodic  boundary 
conditions  (achieved  by  applying  Algorithm  4.1  to  a  larger  interval  with  the  mesh  function  suitably 
extended).  The  numerical  solutions  are  shown  together  with  an  “exact”  solution  computed  on  a 
much  finer  mesh.  These  examples  are  included  to  show  that  the  procedure  can  indeed  handle 
problems  involving  many  interactions. 

The  method  does  have  its  limitations,  of  course.  For  problems  beginning  with  the  initial  con* 
ditions  shown  in  Figure  6.1e,  taking  values  of  kjh  larger  than  10  sometimes  led  to  nonmonotonicity 
of  the  Xy  as  discussed  in  section  5. 
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7.  Several  space  dimensions  and  systems  of  conservation  laws. 

Scalar  conservation  laws  in  several  space  dimensions  can  be  handled  by  means  of  a  spacial 
splitting  or  fractional  step  method.  This  is  a  standard  way  of  extending  one  dimensional  methods 
to  higher  space  dimensions.  Crandall  and  Majda(2]  discuss  the  use  of  fractional  step  methods  for 
computing  weak  solutions  of  conservation  laws. 

In  two  space  dimensions,  the  general  conservation  law  has  the  form 

we(a:,y,0  =  (/Wi,  y,  t))),  +(s(u(i,y,  «)))„■  (7.1) 

Let  5x(A:)  denote  the  approximate  solution  operator  for  the  problem  ui  =  (/(u))x  on  a  time-step  of 
length  k  as  defined  in  the  proceeding  sections.  Similarly,  let  5y(A:)  denote  the  approximate  solution 
operator  for  Ut  =  (^(u))y.  Then  is  computed  by  solving  a  sequence  of  one-dimensional 

problems,  alternating  between  the  x-  and  y-directions.  Using  a  second  order  Strang- type  splitting 
as  described  in  [8],  we  set 

=  Sx(k/2)5y(ifc)5,(ifc/2)u^.  (7.2) 

We  now  have  to  contend  with  errors  both  in  the  ID  difference  scheme  and  in  the  splitting  (7.2). 
ExaLinple  7,1,  Consider  the  problem. 

~  for  0  <  I  <  1,  0  <  y  <  I 

with  initial  data  as  shown  in  Figure  7.1.  This  test  problem  has  been  used  by  Gropp[4j.  The  true 
solution  at  time  1  is  shown  in  Figure  7.2. 

In  this  case  each  ID  problem  involves  a  single  shock,  and  can  be  solved  '‘exactly’’  by  using 
sufficiently  large  time  steps.  The  resulting  shocks  in  the  full  2D  problem  are  still  sharp,  but  are 
ail  parallel  to  the  x-  or  y-axes.  The  result  is  that  shocks  which  have  some  other  orientation,  such 
as  the  "diagonal”  shocks  in  Figure  7.2,  are  represented  by  2ig-zagging  shocks.  The  number  of  zigs 
and  zags  is  directly  proportional  to  the  number  of  ID  problems  solved  to  reach  the  given  time.  In 

this  case  the  best  results  are  obtained  by  using  moderately  small  values  of  t/.  But  we  still  want 

£/  >  I,  or  else  the  ID  procedure  reduces  to  upwind  differencing  and  the  shocks  are  badly  smeared. 
Figure  7.3  shows  some  examples  for  various  values  of  u. 

At  the  present  time  it  is  not  clear  how  to  extend  this  method  to  handle  arbitrary  systems  of 
conservation  laws.  In  general,  a  result  like  Proposition  4.1  will  no  longer  hold,  so  that  interactions 
cannot  be  handled  in  the  same  simple  manner.  For  problems  in  which  the  eigenvectors  of  the 
Jacobian  matrix  df/du  are  nearly  constant  (as  in  some  weak-shock  problems),  it  may  be  possible 
to  handle  interactions  sufficiently  well  to  generalize  the  methods  of  Roe[9].  Research  is  continuing 
on  this  approach. 
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Fij^wro  6.3.  True  and  computed  so/tition*i  Tor  u#  =  iniiini  conditions  from  Figure  6.1b.  fn  fill 

cf»s^:i  h  ~  1/80.  77je  /irnid;LT3  in  pji>‘fr’he:ViJ  infUado  Ow  number  of  Umo^nteps  r.  xded  to  cornpiifo  t.‘ac/i 
so/ution. 


-iVi 


-b- 


Figrirp  Tn/e  iind  rompntcd  sohuhns  for  u,  —  (  -O.lu’).  am/  /nit/a/  conr/itions  from  Fi'^ure  6.1c.  In 
nil  Cftsvs  h  *-=  1/80.  Tiic  nuaiber.i  i/i  pnrcnlheso.s  indicate  the  number  of  time- steps  nvvded  to  eoinpnt ;  each 
solution. 
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Figufv  6.S  Computed  solutions  for  ut  =  (iu^)x  and  initia/  conditions  from  Figure  6, Id  with  h  ^  1/80.  The 
*‘exact"  solution  showii  at  the  bottom  was  computed  using  h  ^  1/240.  The  numbers  in  parentheses  indicate 
the  number  of  time-steps  needed  to  compute  each  solution. 


»c 


Figure  6.6«  Computed  solutions  for  Ut  s:  (— 0.1u^)x  and  initial  conditions  from  Figure  8. id  with  h  ^  1/80. 
The  **exact**  solution  shown  at  the  bottom  was  computed  using  h  1/240.  The  numbers  in  parentheses 
indicate  the  number  of  time-steps  needed  to  compute  each  solution. 


17 


b-ili 


Figure  6.7  Computed  solutions  (or  Ut  ~  initiai  conditions  from  Figure  6.1e  with  h  =  1/80.  The 

"exact”  solution  shown  at  the  bottom  was  computed  using  h  =  1/480.  The  numbers  in  parentheses  indicate 
the  number  of  time-steps  needed  to  compute  each  solution. 


i  ^  0-1 2.S"  i  =  0.2ii 


t  =  o.£ 


Figure  6.8.  Computed  solutions  for  U(  s  ( — 0.1u^)x  and  initial  conditions  from  Figure  G.le  with  h  ^  1/80. 
The  "exact”  solution  shown  at  the  bottom  was  computed  using  h  ==  1/480.  The  numbers  in  parentheses 
indicate  the  number  of  time-steps  needed  to  compute  each  solution. 


II 

u  =  1/2 

U=  -1/2 

u  =  1/4 

u  =  1/2 


u  =  —1/2 


u  =  1/4 


Figure  7.3  Computed  solutions  to  Example  7.1  at  t  —  1.  In  all  cases  h  =  1/25.  The  three  cases  shown  are 
h,k^5h,  and  25h. 
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